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ABSTRACT 

We show three-dimensional magnetohydrodynamical simulations of core col- 
lapse supernova in which the progenitor has magnetic fields inclined to the ro- 
tation axis. The simulations employed a simple empirical equation of state in 
which the pressure of degenerate gas is approximated by piecewise polytropes 
for simplicity. Neither energy loss due to neutrino is taken into account for sim- 
plicity. The simulations start from the stage of dynamical collapse of an iron 
core. The dynamical collapse halts at t = 189 ms by the pressure of high density 
gas and a proto-neutron star (PNS) forms. The evolution of PNS was followed 
about 40 milli-seconds in typical models. When the initial rotation is mildly 
fast and the initial magnetic fields are mildly strong, bipolar jets are launched 
from an upper atmosphere (r ~ 60 km) of the PNS. The jets are accelerated 
to ~ 3 X 10^ km s~^, which is comparable to the escape velocity at the foot 
point. The jets are parallel to the initial rotation axis. Before the launch of the 
jets, magnetic fields are twisted by rotation of the PNS. The twisted magnetic 
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fields form torus-shape multi-layers in which the azimuthal component changes 
alternately. The formation of magnetic multi-layers is due to the initial condi- 
tion in which the magnetic fields are inclined with respect to the rotation axis. 
The energy of the jet depends only weakly on the initial magnetic field assumed. 
When the initial magnetic fields are weaker, the time lag is longer between the 
PNS formation and jet ejection. It is also shown that the time lag is related to 
the Alfven transit time. Although the nearly spherical prompt shock propagates 
outward in our simulations, it is an artifact due to our simplified equation of 
state and neglect of neutrino loss. The morphology of twisted magnetic field and 
associate jet ejection are, however, not affected by the simplification. 

Subject headings: accretion, accretion disks — methods: numerical — MHD — 
supernovae: general 



INTRODUCTION 



Explosion mechanism of core collapse supernova has been an open question for more than 
three decades. Numerical simulations have not succeeded in constructing a convincing model 
of core collapse supernova, although they have been updated and improved steadily. The 
problem is likely to be neither a simple numerical error nor inaccuracy of neutrino transfer. 
It h as been thought that multi-dimensional effects play essential roles in the explosion (see. 



e.g., iBurrows et al.l 120071 . and the references therein). This is based on the fact that the 
spherical symmetric model cannot reproduce explosion while it has been sophisticated to an 
extreme. At the same time, observational evidences have been accumulated for inherent non- 



spherical natures of core c ollapse supernova (see, e.g, IWang et al.ll2002l . l2003l : iHwang et al. 
20041 : [Leonard et al.l 120061 ). 



Magnetic fields and rotation have been thought to be an agent to promote global 
non-sphericity, although i t can be produced w it hout magnetic fields th r ough some other 
mechanisms proposed by iBurrows et al.l (120061 ) . iBlondin &: Mezzacappal (120071 ). and oth- 
ers. As shown by earlier numerical simulations, magnetic fields twisted by rotation pro- 
duces high velocity bipola r jets, if the initial magn etic field is relatively strong and ini- 



tial rotation is fast. Since iLeBlanc fc WilsonI (jl970l ). many magn etohydrodynamica 



ulations of core collapse supernovae have been published (see, e.g. lYamada fc Sawai 



Ardelian et aPbooi 



20061: 



Takiwaki et al. 



Obergaulinger. Aloy &: Miillei 



2004; Sawai 



2006a 



20051: 



Moiseenko et a. 



Obergaulinger et al. 



2006 



a 



20061: 



sim- 



2004 



Shibata et al. 



Burrows et al 



2003). 



However all of them are two-dimensional and have assumed symmetry around the axis. The 
symmetry excludes the possibility that magnetic fields are inclined to the rotation axis. 
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although pulsars are believed to have such magnetic fields. 

In this paper we show three dimensional numerical simulations of core collapse super- 
nova. The initial magnetic field is assumed to be inclined with respect to the rotation axis. 
It is assumed to be stronger than expected from a standard evolutionary model in part 
because a weak magnetic field can have dynamical effects only long afterward and in part 
because it can be strong enough in some circumstances. Such a strong magnetic field may be 
realized in progenitors of magnetars. In a typical model of our simulations, a magnetic torus 
is formed around the PNS and magnetohydrodynamical jets are launched along the initial 
rotation axis. It is also shown that the toroidal component of the magnetic field changes 
its sign alternately in the magnetic torus. We also discuss the dependence on the initial 
magnetic field strength, the initial angular velocity, and initial inclination angle. When the 
initial magnetic fields are weaker, the jets are launched at a later epoch. The total energy 
of the jets depends only weakly on the initial magnetic fields. 

In §2 we summarize our basic model and numerical methods. The results of numerical 
simulations are shown in §3. Discussions are given in §4. Appendix is devoted to the 
numerical schemes that we have developed for the numerical simulations. 



2. MODEL AND METHODS OF COMPUTATION 



2.1. Basic Equations 



As a model of core collapse supernova, we consider gravitational collapse of a massive 
star with taking account of magnetic field. The dynamics is described by the Newtonian 
ideal magnetohydrodynamical (MHD) equations, 

dp 



dt 

dv . ^. 1 
SB 



+ V • {pv) = 0, 



VP 



V xB 



X B 



9^0, 



dt 



= V X (v X B) , 



(1) 
(2) 
(3) 



and 



9 = - V$, (4) 

where p, P, v, B, g, and $ denote the density, pressure, velocity, magnetic field, gravity, and 
gravitational potential, respectively. The gravitational potential, is given by the Poisson 
equation, 

A$ = AnGp. (5) 
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We used the equation of state of iTakahara fc Satol (Il982l ) in which the pressure is ex- 
pressed as 



P = Pc + Pi 



P 



7t - 1 



and 



Pr = 



(6) 
(7) 

(8) 



The index, 7^, is taken to be 1.3. The coefficients, Ki and 7^, are piecewise constant in the 
interval of pj_i < p < Pi- The values are given in Table [TJ The internal energy per unit 
mass is expressed as 

rp p 

e = et + ^dp. (9) 
Jo P 



Accordingly we have the equation of energy conservation, 

d 

— {pE) + V-{pHv) =pv-g, 
where the specific energy (E) and specific enthalpy (H) are expressed as. 



(10) 



V 



Pr 



E = — + -^dp + et + — 



and 



respectively. 



P^ 



P 

H = E + - + — , 
p Stt ' 



Stt 



(11) 



(12) 



2.2. Numerical Grid 

We solved the MHD equations and the Poisson equation simultaneously on a nested grid. 
The nested grid covers a rectangular box of (3.39 x 10'^ km)^ with resolution of Ax = 52.9 km. 
The central eighth volume is covered hierarchically with a finer grid of which cell width 
is the half of the coarse grid. We overlapped rectangular grids of 8 different resolutions 
and achieved very high resolution of Ax = 0.413 km for the central cube of (26.0 km)^. 
The finest grid fully covers the PNS. Since the nested grid has 64^ cells at each level, the 
resolution is roughly proportional to the radius from the center and approximately 4 % 
of the radius. This angular resolution is comparable to that of recent two-dimensional 
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simulations, since the angular resolution is 

2.21 X 10-2 



and A9 



(vr/2)/71 



m 



(7r/2) / 30 = 5.24 x 10"^ in most of them 



Burrows et al.l (120071 ). We call this hierarchically 



arranged grids the nested grid. All the physical quantities are evaluat ed at the cell c enter 
except for the rn a gnetic field. The divergence-free staggered mesh of iBalsaral (1200 ll ) and 
Balsara fc Spicerl (119991 ) is employed for the magnetic field in order to keep V • S = 



wi thin a round-off e r ror. This method is a variant of the constrained transport approach 
of lEvans fc Hawleyl (119981 ). and is optimized for the Godunov-type Riemann solver and 
hierarchical grids. The same type of nested grid has been used for formation of protostars 
from a cloud core (IMatsumoto fc Tomisakal 120041 ). 

The outer boundary condition is set at the sphere of which radius is r = 1.66 x 10^ km. 
The density, pressure, velocity, and magnetic fields are fixed at the initial values outside the 
boundary. 



2.3. Numerical Scheme 



A iRod (jl98ll )-type approximate Riemann solution is employed to solve the MHD equa- 
tions. It takes accouri t of the cold pressure, Pc- Thus it is slightl y diff e rent f rom that of 
Cargo fc Gallicd (119971 ) which is designed to satisfy the property U of lRod (jl98ll ) for an ideal 
gas. The details are given in Appendix. 

We adopted supplementary numerical viscosity to care the carbuncle instability since 
the Roe-type scheme is vulnerable. The supplementa ry viscosity has a larg e value only near 
shocks. The detailed form of the viscosity is given in iHanawa et al.l (120071 ). 



The source term in the equation of energy conservation, pv-g, is evaluated to be the inner 
product of the gravity and the average numerical mass fiux. In other words, we evaluated the 
mass fiux, pv, not at the cell center but on the cell surface. By the virtue of this evaluation, 
the source term vanishes when the mass fiux vanishes. Note that the mass fiux evaluated 
at the cell center may not vanish in a Roe type scheme even when that evaluated on the 
cell surface vanishes. We have found that PNS suffers from serious spurious heating when 
the source term is evaluated from the cell center density and velocity. The spurious heating 
expands PNS to blow off eventually. 



The Poisson equation was solved by the nested grid iteration (IMatsumoto fc Hanawa 



20031 ) as in the simulations of protostar formation. 
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2.4. Initial Model 



Our initial model was constructed from the 15 Mq model of IWoosley et al.l (120021 ) . The 
initial density is increased 10 % artificially to initiate the dynamical collapse. Thus it is 
Po = 6.8 X 10^ g cm~^ at the center. Their model assumes the spherical symmetry and takes 
account of neither rotation nor magnetic field. We have constructed our initial model by 
adding a dipole magnetic field and nearly solid rotation. 



The initial magnetic field is assumed to be 







1 cosO \ 


Be 




-sine 


\ ) 




K ; 



(13) 



in the central core of r < r„ 




fin 



16 



^-^^ cos^ 



16 - — + -| 1 sin^ 







(14) 



in the middle region of < r < 2ra, and 





1 l5Borl 1 


f 2cos0 \ 


Be 


sin 9 




1 8r3 1 


\ B^ ) 




K ) 



(15) 



in the outer region of 2 < r, in the spherical coordinates, where = 846 km. Thus 
the initial magnetic field is uniform inside r < r^, while it is dipolar outside r > Ta- The 
uniform and dipole fields are connected without kink so that the magnetic tension force is 
finite. The electric current density is uniform in the transition region of < r < 2ra in 
this magnetic configuration. 



The initial rotation velocity is expressed to be 

Vq = fi(r) {en x r) 

where 



en 



+ 



- sin 6n 
cos On 



(16) 
(17) 
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and 

/ X \ 

(19) 




in the Cartesian coordinates. The rotation axis is inchned by from the 2;-axis, i.e., from 
the magnetic axis. 

The initial central magnetic field is set in the range of 1.7 x 10^^ G < < 2.0 x 10^^ G 
except in model ROBO. The initial angular velocity is set in the range of 0.31 s^^ < ^^o ^ 
1.21 s~^ except in model ROBO. The models computed are summarized in Table Rl 



The assumed initial magnetic field is stronger than those evaluated by lHeger. Woosley. &: Spruit 



( 120051 ). Our choice is based on the constraint that our three-dimensional numerical simula- 
tions can follow only several tens milliseconds after the bounce. When the initial magnetic 
field is weak, it cannot have any dynamical effects on a relatively short timescale even if 
it is amplified through collapse and rotation. Thus we have assumed rather strong initial 
magnetic field, which can be realized in some progenitors. 



3. RESULTS 



3.1. Non-rotating Non- magnetized Model 



First we show model ROBO having no magnetic field and no rotation as a test of our 
numerical code. In this model the central density increases to reach the maximum value, 
Pmax = 5.71 X 10^^ g cm~^, at t = 187.11 ms. It settles down to the equilibrium value, 
Peq = 4.66 X lO^'' g cm~^, after several times of radial oscillation. The oscillation period is 
1.25 ms. Radial shock waves are initiated by this core bounce. The first one, the prompt 
shock wave, reaches r = 595 km at t = 222.71 ms, where r denotes the radial distance from 
the center. It reaches the boundary of computation (r = 1700 km) at t = 303 ms, while 
the expansion velocity decreases. The fast propagation of prompt shock is an artifact due 
to our simplified EOS. If we had taken neutrino transfer into account, the prompt shock 
should have stalled aro und r ^ 100 kra roughly within 100 ms after the bounce (see, e.g. 



Sumiyoshi et al.ll2005l : iBuras et al.l 120071 : iBurrows et al.ll2006l ). 



Figure [T] shows the radial density profiles at t = 189.3, 210.9, and 232.2 ms in model 
ROBO. The density decreases monotonically with increase in the radius. The density gradient 
is steep around r ~ 20 km. In the outer region of r > 30 km, the density decreases gently 
and roughly proportional to r~^. Thus we regard the layer of p = 10^^ g cm~^ as the surface 
of PNS for simplicity. The mass and radius of PNS are 1.10 Mq and 26.6 km, respectively. 
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at t = 210.9 ms. The mass increases by 4.1x10 ^ M© between t = 189.3 and 232.2 ms. 

Non-radial oscillation has not been excited, although our numerical code does not as- 
sume any symmetry. This is most likely to be due to short computation time, i.e., only 
40 ms after th e bounce. Note that t he i = 1 mode becomes appreciable around 200 ms after 
the bounce in lBurrows et al.l (120061 ). 



3.2. Typical Model with Inclined Magnetic Field 

In this subsection we show model R12B12X60 as a typical example. The initial rotation 
period is small compared to the free-fall timescale, fio/V^TrGpo = 1-59 x 10"^, and initial 
rotation energy is much smaller than the gravitational energy (|-Ekin/-E'grav| = 5.0 x 10^^). 
The magnetic energy is also much smaller than the gravitational energy (|-E'mag/-E'grav| = 
2.9 X 10"^). 

Figure [2] shows the evolution of central density (pc) for the period of 180 ms < t < 
230 ms. The central density reaches its maximum, 5.49 x 10^'* g cm~^, at t = 189.05 ms as 
well as in model ROBO. The period of dynamical collapse is a little longer and the maximum 
density is a little lower. The rotation and magnetic field delays the collapse a little as has 
been shown in earlier simulations. The PNS is only slightly flattened by rotation. 

At t = 190.04 ms (slightly after the bounce), the PNS has angular velocity of Qc = 
6.02 X 10^ s~^ and magnetic field of Be = 7.56 x 10^^ G at the center. The angular velocity 
and magnetic field increase by a factor of 5.0 x 10^ and 3.8 x 10^, respectively, from the 
initial values, while the density increases by a factor of 6.1 x 10^. These enhancements 
are consistent with the conservation of the specific angular momentum and magnetic flux, 
since the collapse is almost spherical. The angular velocity and magnetic field increase in 
proportion to p^^^ when the collapse is spherical. The rotation axis changes little. At this 
stage the centrifugal force is only 3 % of the gravity at the center of the PNS. The magnetic 
force is much weaker than the gravity and than the centrifugal force. 

The change in the magnetic field is illustrated in Figure [31 The magnetic field is almost 
radial near the end of dynamical collapse as shown in the top panels, in which the purple 
lines denote the magnetic field lines at t = 188.28 ms by bird's eye view. This is because 
the magnetic field is stretched in the radial direction by the dynamical collapse. The radial 
component of the magnetic field, B,., is positive in the upper half of z > 0, while it is 
negative in the lower half. Thus the split monopole is a good approximation to the magnetic 
field at this stage. 



- 9 - 



The magnetic field is twisted by the spin of PNS as shown in the middle and bottom 
panels of Figure [31 The azimuthal component of the magnetic field is amplified to have a 
large amplitude in the upper atmosphere of the PNS (9 km ^ t < 14 km). The increase in 
the azimuthal component decreases the plasma beta down to (3 = -Pgas/-Pmag — 0.03. The 
azimuthal component of the magnetic field is small inside the PNS, since the angular velocity 
is nearly constant. It is also small in the region very far from the center (r > 60 km) since 
the angular velocity is very small. It is also small near the rotation axis. Thus the region of 
strong twisted magnetic field has a torus-shape. 

The structure of twisted magnetic field is different from that in an aligned rotator. 
When the initial magnetic field is aligned to initial rotation axis, the twisted magnetic field 
has opposite directions in the upper and lower halves. The azimuthal component vanishes on 
the equator of rotation and has a large amplitude in the upper and lower tori. The magnetic 
field is uni-directional in each torus. In case of oblique rotator, however, these tori are mixed 
into a torus, in which the azimuthal component of magnetic field is bi-directional as a result. 
Figure m shows the distribution of toroidal component of magnetic field, = (e<^ ■ B), at 
t = 197.92 ms, where = bq x r/\e^ x r\ [cq denotes the unit vector along the initial 
rotation axis and heads upper left in Figure H]). The toroidal component changes its sign 
alternately with an average interval of ~ 5 km. 

Figure [5] is the same Figure H] but for the 4.89 ms later stage. The magnetic field is 
wound more tightly inside the PNS and the toroidal component has extended outward. 

A similar magne tic field is obtained semi-analytically as a model of pulsar magnetosphere 



by iBogovalovi (119991 ). He approximated the initial magnetic field by split-monopole one and 
considered oblique rotation. As shown in his Figure 4, the toroidal component changes its 
sign with a regular interval in his pulsar wind solution. His idealized magnetic configuration 
is realized in our simulation. The magnetic multi-layers is inevitably formed when the initial 
magnetic field is split-monopole-like and inclined with respect to the rotation axis. 

Figure [6] shows later evolution of the magnetic field. The torus of twisted magnetic field 
expands slowly and bipolar jets are launched along the rotation axis. The jet reaches r = 
400 km at the last stage of computation [t = 228.99 ms). The jet velocit y exceeds 3 x 10^ km 



s ^. T he jets are driven mainly by the magneto-centrifugal mechanism of iBlandford fc Payne 



(119821 ). Figure [7] shows the evolution of rotation velocity around the initial rotation axis. The 
rotation is nearly rigid in the sphere of r < 10 km just before the bounce {t = 188.52 ms). 
The rotation velocity increases up to ~ 5 x 10^ km s^^ by the magnetic torque near the 
rotation axis. The magneto-centrifugal force is strong enough to drive jets. Although the 
centrifugal force is perpendicular to the rotation axis, the component perpendicular to the 
magnetic field is cancelled by the strong magnetic force. Thus the gas is accelerated along 
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the poloidal magnetic field, i.e., along the rotation axis by the Blandford-Payne mechanism. 

The increase in the radial velocity follows that in the rotation velocity. The acceleration 
of jets are shown in Figure [HI The radial velocity is still low at the stage shown in the upper 
left panel (t = 207.56 ms). The high velocity jets emerge not from the PNS surface but from 
the upper layer of r ~ 60 km. The mass flux through the sphere of r = 300 km is M = 0.0, 
7.64, 5.43, 5.58, 3.75, 2.07 Mq s"^ at t = 190.75, 207.56, 212.96, 222.80, and 228.99 ms, 
respectively. Figure [H] shows the jets by bird's eye view at the stage of t = 222.80 ms. The 
jets are bipolar and well coUimated. 

The magnetic force is comparable to the gas pressure in the jets. Figure [TO] denotes the 
magnetic pressure distribution at t = 208.42 ms. The thick solid curve denotes the magnetic 
pressure, \B\'^/87i, along the initial rotation axis, while thin solid curve does that on the 
equator. The magnetic pressure is enhanced by winding in the range 50 km ^ r < 100 km 
on the axis. Also the dynamical pressure (rotation energy) is enhanced in the same region. 
Since these energies are large enough, the jets will extend outwards even if the prompt shock 
has stalled around r ~ 100 km. 

We computed the energy of magnetic field, rotation and jets to evaluate the efficiency 
of energy conversion. The magnetic energy distribution is evaluated by 

emag(r, t) = J |r3M^l_^Ll^ sin^d^d^, (20) 

i.e., the energy stored in a unit logarithmic radial distance. The total magnetic energy is 
expressed as 

^mag(i) = j £mag(r, t) dlur . (21) 

Figure [TT] shows that the magnetic energy has a sharp peak of ^mag = 1-88 x 10^^ erg in the 
layer of r ~ 18 km at t ~ 196 ms. The peak of the magnetic energy shifts outward at the 
apparent radial velocity of 3 x 10^ km s~^, which coincides with the Alfven velocity. The 
peak declines beyond r > 80 km. 

Figure [12] is the same as Figure [11] but for the radial kinetic energy stored in a unit 
logarithmic radial distance, 

£kin,rad(r, t) = r 81116(16 dip, (22) 

where n = r/\r\ denotes the unit radial vector. In the region far from the center, the dy- 
namical collapse dominates in the radial kinetic energy. After the bounce {t > 189.72 ms), 
the prompt shock propagates and the region of the dynamical collapse retreats. (The prop- 
agation of prompt shock is mainly due to neglect of neutrino loss. If we had incorporated 
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the neutrino cooling, the prompt shock should have stalled around 100-200 km.) The radial 
kinetic energy is small in the region of r < 60 km after the bounce. The bipolar jets increase 
the radial kinetic energy in the region of r > 80 km. The rise in ^kin.rad coincides with the 
decline in e^^g- This is an evidence that the jets are accelerated by magnetic force. 

Figure [13] shows the distribution of the kinetic energy stored in a unit logarithmic radial 
distance, 

ekm,rot(r, t) = / r sm ededip. (23) 

A large amount of rotation energy is stored in the region of 10 km ^ t < 20 km. Only a 
small fraction of it is converted into the energy of twisted magnetic field and eventually into 
the energy of jets. 

Figure [14] shows the evolution of energy stored in the volume of r < 63 km for each 
component. The gravitational energy, which is evaluated to be 

is the most dominant and the internal energy is comparable. The thick solid curve denotes, 
Aii^grav = -Egrav — -Sgrav,min, the difference from the minimum value, i.e., the gravitational 
energy at the stage of the maximum central density (bounce). The rotation energy is order 
of magnitude smaller than them and the magnetic energy is further smaller. Only a small 
fraction of the rotation energy is converted into magnetic energy, which is mostly due to 
toroidal magnetic field. The energy available for jet ejection is limited by the conversion 
factor from rotation energy to magnetic energy. The radial kinetic energy is only of the 
order of ~ 10^^ erg in the period t > 195 ms since the prompt shock and jets are outside 
the region. For comparison we show the evolution of the energy stored in model ROBO by 
thin curves. 



Note that the rotational energy available is much smaller than those in lObergaulinger et al 



(l2006al ). Since the initial rotation energy was much larger in their simulation, the PNS 
shrunk appreciably after liberating the angular momentum through magnetic braking. Even 
if the rotation energy were released completely in our model, the PNS would not shrink 
appreciably. 

When the fast jets are ejected along the initial rotation axis, two other radial flows are 
observed. Figure [15] shows the velocity distribution on the plane of x = at t = 229.77 ms. 
One is slow outflow extending near the equator of initial rotation. The outflow velocity is 
approximately 2.5 x 10^ km s"^. The other is fast radial inflow located between the jets and 
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equatorial outflow. The inflow is less dense and its dynamical pressure is much smaller than 
the magnetic pressure. 

As shown earlier, the bipolar jets emanate from r ^ 60 km. In the outer region of 
r > 60 km, the density is lower than p < lO^'' gm cm~^ (see Figure [1] for the average 
density distribution in model ROBO) and hence the Alfven velocity is high. In other words, 
the magnetic force dominates over the pressure force. The centrifugal force is also important 
in the outer region. If the magnetic field corotates with the PNS, the rotation velocity is 
evaluated to be 



3.6 X 10^ 



km s ^ , 



6 X 103 s-i J V60 km/ 

where w denotes the distance from the rotation axis. The rotation velocity is close to the 
Kepler ian velocity. 



GM 



Vk 



4.70 X 10^ 



M 



1/2 



60 km 



-1/2 



km s 



(26) 
(27) 



where M and r denote the PNS mass and the distance from the center, respectively. In other 
words, the centrifugal force is comparable with the gravity. Thus the foot point of MHD jets 
coincides with the inner edge of the region in which the magnetic and centrifugal forces are 
dominant over the gravity. 

As shown in Figure H] the twisted magnetic field are ordered and has no structures sug - 
gesting development of magneto-rotational instability (MRI; see e.g.. lAkiyama et al.ll2003l ). 
Howe ver, this is likely due to li r aited spatial resolution and will not exclude the possibility of 
MRI. lEtienne. Liu. &: Shapiro! (120061 ) demonstrated that MRI can not grow unless the cell 
width is shorter than a tenth wavelength of the fastest growing mode. Since the wavelength 
is evaluated to be 



A 



MRI 



47r 



1.18 



^1/2 



B 



dr J ^/4vrp 
d^^/d\nr T^/^ 
(3000 s-i)^ 



B 
1015 G 



10 



14 



P_ 

g cm" 



-1/2 



km , 



(28) 
(29) 



we need the spatial resolution of ~ 120 m to observe MRI. 
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3.3. Dependence on Bq 

To examine the effect of initial magnetic field we made 6 models by changing only Bq 
from model R12B12X60. Figure [16] shows the maximum radial velocity (fr,max) as a function 
of time. It declines sharply in the period of t < 195 ms, since the prompt shock wave slows 
down. The early decline is delayed when Bq is larger. The delay is, however, smaller than 
1 ms since the magnetic energy is much smaller than the gravitational energy at the PNS 
formation. See Table [X] for comparison of models at the bounce. 

The late rise in Wr,max, i-e., the launch of MHD jets depends strongly on Bq. When Bq 
is larger, the radial velocity rises earlier and stays at a high level. When Bq is smaller than 
1.0 X 10^^ G, the radial velocity increases late but decreases soon before reaching a high level. 

When Bq is large, the magnetic field is less tightly twisted since the twisted component 
drifts upward faster. Figured?] is the same as Figure H] but for model R12B16X60, in which 
Bq is 4/3 times larger than in the standard model R12B12X60. The toroidal component, 
B^, changes its sign with a longer average interval of 5-6 km on the average, while it does 
with an interval of 5 km in the standard model. At t ~ 206 ms, the twisted magnetic field 
reaches r = 60 km as shown in Figure [18] and the MHD jets initiate. 

On the other hand, the magnetic field is more tightly twisted when Bq is smaller. 
Figure [19] is the same as Figured] but for model R12B8X60. The toroidal component, B^p, 
changes its sign with a shorter interval of 3-4 km on the average. The twisted magnetic field 
rises up slowly and dwindles as shown in Figure [20] Accordingly the jets are launched but 
late and weak as shown in Figure [THl where the evolution of the maximum radial velocity, 
Vr,ma.x, IS showu for various models having different Bq. The weakness of the jet is at least 
partly due to numerical diffusion. Remember that the spatial resolution is 0.826 km in 
the central cube of (53.0 km)^. Thus the numerical diffusion is appreciably large for the 
magnetic multi-layers since the typical interval is less than 10 km. The MHD jets would be 
more powerful if the numerical diffusion were suppressed. 

See Table [A] to compare models at the final stages. 

3.4. Dependence on ^Iq 

To examine the effect of initial rotation, we made 5 models by changing only Qq from 
model R12B12X60. All the models show qualitatively similar results. The differences are 
mainly quantitative. 

The initial rotation is twice slower in model R6B12X60 than in model R12B12X60. 
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Accordingly the PNS has twice lower angular velocity in model R6B12X60. The magnetic 
field is twisted by rotation also in model R6B12X60 but the toroidal component is weaker 
since the rotation is slower. The MHD jets are launched but at a little later epoch. 

Figure [21] shows the evolution of the maximum radial velocity, fr,max, as a function 
of time for various models having different Qo- The early decline in fr,max is due to the 
deceleration of the prompt shock. The late rise in Vr^max is due to the launch of jets. When 
Qq is larger, the PNS is formed a little later and the jets are ejected a little earlier. The 
maximum radial velocity is slower when Qq is small. 

The rotation energy of PNS is proportional to Q^. It is i?rot = 3.8 x 10^^ erg in model 
R18B12X60 while it is E.^t = 4.5 x 10^^ erg in model R6B12X60. The energy of the jets, 
which is evaluated to be radial kinetic energy of outflowing gas, is also large in a model 
having a large ^o- 

3.5. Dependence on 6q 

To examine the effect of the initial inclination angle, we made 5 models by changing 
only 6^ from model R12B12X60. Also in these models, the MHD jets are launched along 
the initial rotation axis (see Fig. [2^ . 

Figure [23] denotes the magnetic multi layer formed in model R12B12X30. The structure 
of magnetic multi layer depends on the inclination angle. When the inclination angle is 
smaller, it is conflned in a narrower region around the equator. The radial interval of 
changing depends little on the inclination angle. 

Figure [2l] shows the evolution of the maximum radial velocity in these models. When 
On is larger, the maximum radial velocity rises earlier. In other words, the MHD jets are 
launched earlier when the rotation axis is inclined. 

Although the rise is different, the maximum radial velocity reaches a certain value 
independently of the inclination angle. In other words, the terminal velocity is independent 
of the inclination angle, 6q. 

4. SUMMARY AND DISCUSSIONS 

In this paper we have shown three dimensional MHD simulations of core collapse super- 
nova for the flrst time. The numerical simulations explore the effects of inclined magnetic 
fleld and dependence on the initial magnetic fleld and rotation. We summarize the results 
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and discuss the implications. 

First we have confirmed that the MHD jets are ejected along the initial rotation axis. 
This is because the energy of rotation dominates over the magnetic energy at the moment 
of PNS formation. The magnetic force is too weak to change the rotation axis appreciably. 
Thus the magneto-centrifugal force accelerates gas along the initial rotation axis through 
the Blandford-Payne mechanism (see §3.2). 

If the initial magnetic field were relatively strong, the rotation axis could change appre- 
ciably and also the je t direction could be different. A similar problem has been studied by 
Machida et al. for formation of a protostar. They considered collapse of a rotating 



molecular cloud core having an oblique magnetic field. They have found that the evolution of 
magnetic field and rotation axis depends on their relative strength. When the magnetic field 
is relatively strong, the magnetic braking acts to align the rotation axis with the magnetic 
field. T hen the jets are ejected in the d irection parallel to the initial magnetic field as shown 



first by lMatsumoto fc Tomisakal (120041 ). 



The relative strength between rotation and magnetic field can be evaluated from the 
ratio of the angular velocity to the magnetic field, VL/ B. The ratio remains nearly constant 
during the dynamical collapse since the free-fall timescale is very short. Both the angular 
velocity and magnetic field increase in proportion to the inverse square of core radius, since 
both the specific ang ular momentum and m agnetic flux change little during the short dynam- 



ical collapse phase. iMachida et al.l (120061 ) proposed VLq/Bq > 0.39 G^'"^ , clS 8b criterion 
for the jets parallel to the initial rotation axis, where Cg denotes the isothermal sound speed 
of the molecular cloud. Since the dynamics of collapse is similar, we can expect that the 
criterion holds also for core collapse supernova if we replace Cg with an appropriate one. The 
criterion is rewritten as 



iln c 







GV^ Bo 

The left hand side is evaluated to be 



> 0.39. (30) 



3.87r^U^) i,^,^,]. (31) 



G^/^Bo • \ls-^J\W^Gj VlO^kms 

The criterion is consistent with our numerical simulations since the sound speed increases 
from Cg = 10^ km s"'^ to 10^ km s"'^ during the dynamical collapse. Although the assumed 
initial magnetic field is strong, it is still too weak to change the rotation axis unless the initial 
rotation is slow. Thus it is reasonable that the rotation is unchanged during the dynamical 
collapse since young pulsars are spinning fast. 

Next we discuss the fate of magnetic multi-layers in which the toroidal magnetic field 
changes its direction with a regular interval. The magnetic multi-layers are a natural outcome 
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of oblique rotation as shown in the previous section. These layers are potentially unstable 
against reconnection, although no features are seen for reconnection. It is also interesting to 
study the magnetic multi-layers with a higher spatial resolution. If the magnetic fields are 
reconnected, a large amount of the magnetic energy is released to lead an explosive process. 

Next we discuss the lag between the bounce an d jet ejection. When the initial magnetic 



f ield i s weak er, the lag is lo n ger a s already shown in lArdeljan et al.l (120041 ) , iMoiseenko et al. 



( 120061 ) , and iBurrows et al.l (120071 ) . Our numerical simulations have suggested that the lag 
is related to the Alfven transit time; the MHD jets are ejected when the twisted magnetic 
field reaches a certain radius, i.e., 60 km in our simulations. When the initial magnetic field 
is weak, the Alfven transit time is longer and the magnetic field is amplified for a longer 
duration before the jet ejection. Since the larger amplification compensates for a weak seed 
field, almost the same amount of toroidal magnetic field is generated irrespectively of the 
initial magnetic field strength. This implies that strong MHD jets can be ejected even if 
the initial magnetic field is weak. If we could suppress numerical diffusion by improving 
resolution, strong MHD jets should be ejected also in model R12B5X60 and others having a 
weaker initial magnetic field as discussed in the previous section. 

The Alfven transit time is evaluated to be 

TA = / —dr (32) 
J Va 



—— dr , (33) 

where rj (~ 60 km) denotes the radius at the foot point of the jets. The radial magnetic 
field decreases in proportion to the inverse square of the radius (oc r^^) since the split 
monopole is a good approximation for the magnetic field at the epoch of PNS formation. 
The density decreases also roughly proportional to the inverse square of the radius (oc r~^) in 
the upper atmosphere of the PNS (see Figured]). Accordingly the Alfven velocity is roughly 
proportional to the radius, va oc r and the Alfven transit time depends only weakly on Vj. 

The above argument suggests a factor controlling the jet energy. The rotation energy of 
the PNS is much larger than the energy of jets. If the magnetic field is twisted for a longer 
period, i.e., if the Alfven transit time is longer, a larger fraction of the rotation energy is 
converted into the jet energy. The Alfven transit time can be extended if the magnetic field 
is twisted in a deep interior of the PNS. 
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A. 



An Approximate Riemann Solution of the MHD Equations for Non-ideal 

Equation of State 



For the approximate EOS of iTakahara &: Satol (Il982l ). we have found a numerical flux 



which satisfies t h e pro perty U of Roel ( 1981 ). The numerical flux is the same as that of 



Cargo fc Gallicd (119971 ) except for the correction term. 



G 



where 



^cold + P ■ 



^cold 



cold 



It 



1 Ap 



P, 



cold 



P^ 



dp . 



(Al) 



(A2) 



Here the bared symbols denote the Roe average while the symbols with the capital delta de- 
note the differences between the two adjacent cells. The average specific enthalpy, should 
be replaced with H — G in the computation of a and a^,. Accordingly, the characteristic 
speeds propagating in the x-direction are evaluated to be 



{it 



H - G 



Anp 



6b^ 



(A3) 



where 



(7* - 1; 

-ht-2 



H - G 



Bl + Bl + Bl 



It - 1 



AB'y + ABl 



6b' 



(A4) 



(A5) 



7t- 2 47r(ypr+ypi)2- 

Here, {u, v, w) and {B^, By, Bz) denote the velocity and magnetic field, respectively. The 
symbols, ABy and Ai?^, denote the difference between the two adjacent cells, while the 
symbols, pl and pr denote the density in the left hand side and that in the right hand side, 
respectively. 
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At the same time, the correction term G should be added to the last component of 
the eigen vector of the en tropy wave. This correction term is similar to that obtained by 



Nobuta fc Hanawal (119991 ) for the numerical flux of hydrodynamical equations. Thus the 



right eigen vector for the entropy wave is expressed as, 

t I . - - - ^ ^ ^ + v'^ + w"^ 



^entropy ( 1; ^) 

0, 0, 0, ^ + G] . (A6) 
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Fig. 1. — The angle averaged density, p (r, t), is shown as a function of r for model ROBO. 




m w- ^ 



Fig. 2. — The central density is shown as a function of time for model R12B 12X60. 
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Fig. 3. — The evolution of magnetic field in model R12B12X60. Each panel denotes the 
magnetic field hnes (purple lines) and the isovelocity surface of radial velocity, f^, at a given 
stage. The panels are arranged in the time sequence from top to bottom. The top panels 
denote the stage of t = 188.28 ms, while the bottom ones do that of t = 191.10 ms. The 
left panels denote the central cube of (1691 km)^, while the central and right ones do the 
zoom-up views of of (423 km)^ and (26 km)^, respectively. The left color bars denote the 
radial velocity on the surface of the cubes in the unit of 10'^ km s~^. 
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Fig. 4. — Structure of magnetic torus is shown by three cross sections and a bird's eye view. 
The color denotes the azimuthal component of magnetic field, B^p = e^p ■ B on the planes 
of x = (upper left), y = (upper right), and z = (lower right) at t = 197.92 ms. The 
contours denote in unit of 10^^ G. The lower right panels magnetic field lines (purple) 
and isosurface of B^p by the bird's eye view. 




Fig. 5. — The same as Fig. Hlbut for the stage at t =202.81 ms. 




Fig. 6. — The same as Fig. [3] but for later stages of t = 192.90 ms, 193.97 ms, 208.70 ms, 
and 228.99 ms. 
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Fig. 7. — The evolution of rotation velocity. Each panel denotes the distribution of 
log v^p (km s~^) around the initial rotation axis in the cross section x = in model 
R12B12X60 by color and contours. The arrows denote the velocity within the plane. They 
denote the stages at t = 188.52, 190.75, 207.56, and 212.96 ms. 
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Fig. 8. — The radial velocity [vj.) distribution is denoted by color and contours for the stages 
at t = 207.56, 212.57, 222.80, and 228.99 ms in model R12B12X60. The numbers attached 
to the contours denote Vr in unit of 10^ km s~^. The arrows denote the velocity within the 
cross section. 
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Fig. 9. — Jets and fast radial inflows in model R12B12X60. The panel denotes the central 
cube of (423 km)^ at t = 222.80 ms by bird's eye view. The magnetic field lines are denoted 
by the purple lines, while the radial velocity is denoted by the isosurfaces. The blue denotes 
the fast radial inflow {vr < —1.0 x 10^ km s^^), while the others do jets {v,- > 2.0 x 10^ km 
s~^). The color bar is for the radial velocity distribution on the cube surfaces. 




Fig. 10. — Pressure distribution is shown for the stage at t = 208.42 ms. The solid curves 
denote the magnetic pressure, \B\'^/8n, while the dashed curves denote the gas pressure. 
The dash-dotted curves denote the dynamical pressure, p|i^p. The thick curves are for the 
values on the initial rotation axis while the thin curves for those on equator. 
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Fig. 11. — The magnetic energy stored in a unit logarithmic radial distance, e^^gif^, t), is 
shown by darkness. The contour levels are set to be A^mag = 10''^ erg. The abscissa is the 
time (t) while the ordinate is radial distance (r). 
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Fig. 12. — The same as Fig. [TT]but for the radial kinetic energy, £:kin,rad- 
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Fig. 13. — The same as Fig. [TT]but for the rotational energy, ekin,rot- 
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Fig. 14. — The evolution of energy in model R12B12X60. Each thick curve denotes a 
component of the energy stored in the sphere of r < 63 km. The thin curves denote those 
in model ROBO. 



Fig. 15. — The velocity distribution are denoted by the arrows on the plane of a; = at 
t = 229.77 ms. The scale is shown on the upper left corner. The darkness denotes the 
region of positive Vr while the contours denote isodensity curves and labeled by log p. 




Fig. 16. — The maximum radial velocity is shown as a function of time for various models 
having different initial magnetic field. 




Fig. 17.— The same as Fig. 1 but for model R12B16X60. 
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Fig. 18.— The same as Fig. [ID but for model R12B16X60. 




Fig. 19.— The same as Fig. 1 but for model R12B8X60. 
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Fig. 20.— The same as Fig. [ID but for model R12B8X60. 
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Fig. 21. — The same as Fig. [16] but various models having different initial angular velocity. 



-37- 



J.. ■ ■«■•■" 
J. 

■I r 




I °^ • 

V 




I 



--."..■aV. ■ 



Fig. 22. — Comparison of radial velocity distribution between the models having different 
rotation axis (On)- The inclination angle is 6^ = 15° (top), 30° (middle), and 60° (bottom), 
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Fig. 23.— The same as Fig. 1 but for model R12B12X30 at t =202.79 ms. 
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Fig. 24. — The same as Fig. EHfor various models having different initial inclination angle. 
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Table 1: Model parameters for P^. See Eq. ID 



i Pi Ki 7i_ 

1 4.0 X 10^ 1.767 X 1027 4/3 

2 1.0 X 10^2 2.446 X lO^o 1.31 

3 2.8 X 10^^ 4.481 x 10^3 4/3 

4 1.0 X 10^5 1.080 X 10^5 2.5 
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Table 2. Summary of the models. The initial gravitational energy is 
W = —4.75 X 10^^ erg in all the models. The models underlined are shown in detail in the 

main text. 



model 


^^0 


-E'kin 


|-£'kin/-E'grav| 


Bo 


-S'mag 


|-^mag 




;rav | 


en 






(erg) 


(%) 






(G) 


(erg) 


(%) 








ROBO 


































R12B12X60 


1.2 


2.2x10^^ 


5.0 


X 


10- 


-4 


2.0x10^2 


1.4x10^^ 


2.9 


X 


10- 


-4 


60° 


R18B12X60 


1.8 


4.9x10^8 


1.1 


X 


10- 


-3 


2.0x10^2 


1.4x10^^ 


2.9 


X 


10- 


-4 


60° 


R15B12X60 


1.8 


3.4x10^8 


1.1 


X 


10- 


-3 


2.0x10^2 


1.4x10^8 


2.9 


X 


10" 


-4 


60° 


R8B12X60 


0.81 


9.6x10^^ 


2.0 


X 


10- 


-4 


2.0x10^2 


1.4x10^^ 


2.9 


X 


10- 


-4 


60° 


R6B12X60 


0.61 


5.4x10^^ 


1.1 


X 


10" 


-5 


2.0x10^2 


1.4x10^^ 


2.9 


X 


10- 


-4 


60° 


R3B12X60 


0.31 


1.4x10^^ 


2.0 


X 


10- 


-4 


2.0x10^2 


1.4x10^^ 


2.9 


X 


10- 


-4 


60° 


R12B16X60 


1.2 


2.2x10^^ 


5.0 


X 


10- 


-4 


2.7x10^2 


2.5x10^^ 


5.8 


X 


10- 


-4 


60° 


R12B8X60 


1.2 


2.2x10^^ 


5.0 


X 


10- 


-4 


1.3x10^2 


6.2x10^^ 


1.4 


X 


10- 


-4 


60° 


R12B6X60 


1.2 


2.2x10^^ 


5.0 


X 


10- 


-4 


1.0x10^2 


3.5x10^^ 


8.1 


X 


10- 


-5 


60° 


R12B5X60 


1.2 


2.2x10-^^ 


5.0 


X 


10- 


-4 


8.4x10" 


2.4x10^^ 


5.6 


X 


10- 


-5 


60° 


R12B4X60 


1.2 


2.2x10^^ 


5.0 


X 


10- 


-4 


6.7x10" 


1.6x10^^ 


3.6 


X 


10- 


-5 


60° 


R12B1X60 


1.2 


2.2x10^8 


5.0 


X 


10- 


-4 


1.7x10" 


9.7x10^^ 


2.3 


X 


10- 


-6 


60° 


R12B12X30 


1.2 


2.2x10^8 


5.0 


X 


10- 


-4 


2.0x10^2 


1.4x10^^ 


2.9 


X 


10- 


-4 


30° 


R12B12X15 


1.2 


2.2x10^^ 


5.0 


X 


10- 


-4 


2.0x10^2 


1.4x10^^ 


2.9 


X 


10- 


-4 


15° 


R12B12X7 


1.2 


2.2x10^^ 


5.0 


X 


10- 


-4 


2.0x10^2 


1.4x10^^ 


2.9 


X 


10- 


-4 


7° 


R12B12X0 


1.2 


2.2x10^^ 


5.0 


X 


10- 


-4 


2.0x10^2 


1.4x10^^ 


2.9 


X 


10- 


-4 


0° 
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Table 3. Comparison of models at the bounce. 



J 1 

model 


t 




T/\W\ 


M/\W\ 


I 


M 


^max 


-'-'max 




(ms 


) 


(%) 


(%) 


[10°" erg) 


(10^° erg) 


/I TT \ 

(kHz) 


[10 ° G) 


ROBO 


191, 


.9 


0, 


.5854 


0. 


,000 


8.359 


0, 


.000 


0.1421 


0.000 


R12B12X60 


192, 


.2 


2, 


.137 


0. 


.1002 


30.46 


1 


.428 


29.97 


17.96 


R18B12X60 


192, 


.2 


4, 


.315 


0. 


.07094 


57.21 


0, 


.9405 


38.21 


13.48 


R15B12X60 


194, 


.1 


2, 


.703 


0. 


.1679 


37.72 


2, 


.344 


38.44 


21.25 


R8B 12X60 


191 


.9 


1 


.309 


0. 


.07185 


18.61 


1 


.021 


26.50 


16.08 


R6B12X60 


193, 


.1 


0, 


.9306 


0. 


,07122 


13.24 


1, 


.013 


16.72 


18.46 


R3B 12X60 


193, 


.2 


0, 


.6340 


0. 


,04088 


9.019 


0, 


.5815 


6.329 


18.04 


R12B16X60 


192, 


.7 


2, 


.033 


0. 


,1721 


28.91 


2, 


.447 


28.20 


23.09 


R12B8X60 


192, 


.3 


2, 


.150 


0. 


,05520 


30.53 


0, 


.7840 


28.39 


13.41 


R12B6X60 


192, 


.6 


2, 


.022 


0. 


.1000 


28.17 


1 


.394 


7.656 


16.68 


R12B5X60 


192, 


.3 


2, 


.180 


0. 


.02334 


30.96 


0, 


.3314 


27.82 


9.320 


R12B4X60 


220, 


.7 


2, 


.887 


0. 


.01407 


41.06 


0, 


.2001 


31.03 


9.416 


R12B1X60 


192, 


.5 


2, 


.127 


0. 


.08790 


30.23 


1 


.249 


35.31 


16.50 


R12B12X30 


194, 


.9 


1, 


.986 


0. 


,1384 


27.98 


1, 


.950 


42.41 


16.26 


R12B12X15 


194, 


.9 


1, 


.981 


0. 


,1423 


27.91 


2, 


.005 


41.17 


16.83 


R12B12X7 


192, 


.7 


2, 


.036 


0. 


.03528 


28.22 


0, 


.4889 


43.35 


6.467 
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Table 4. Comparison of models at the final stage. 



J 1 

model 


t 

(ms 


) 


Tl\W\ 
(%) 


MI\W\ 
(%) 


T M 

(10^° erg) (10^° erg) 


(10^° erg) (kHz) 


-^max 

(10^5 G) 


ROBO 


232. 


,2 


0, 


.1552 


0. 


,000 


2.312 


0, 


,000 


61.53 


5.123 


0.000 


R12B12X60 


229. 


,0 


2, 


.277 


0. 


,1835 


32.44 


2, 


.615 


1.002 


27.80 


12. 


,80 


R18B12X60 


215. 


,8 


3 


.497 


0. 


,1955 


48.38 


2, 


.705 


57.54 


32.85 


15. 


,89 


R15B12X60 


224. 


,4 


2, 


.988 


0. 


,2084 


42.11 


2, 


.937 


13.59 


32.59 


14. 


,21 


R8B 12X60 


228. 


,7 


1 


.200 


0. 


,1547 


17.33 


2, 


.235 


23.50 


19.85 


13. 


,99 


R6B12X60 


231, 


,9 


0, 


.7440 


0. 


,1066 


10.94 


1, 


,568 


47.14 


12.32 


14. 


,89 


R3B12X60 


233. 


,2 


0, 


.2592 


0, 


,06470 


3.853 


0, 


,9616 


63.81 


5.220 


11. 


,56 


R12B16X60 


230. 


,5 


2, 


.010 


0. 


,2382 


28.77 


3, 


,409 


8.921 


24.99 


12. 


,74 


R12B8X60 


227. 


,7 


2, 


.514 


0. 


,1548 


35.67 


2, 


.197 


1.261 


38.35 


13. 


,37 


R12B6X60 


224. 


,9 


1 


.773 


0. 


,1817 


23.26 


2, 


.383 


81.68 


12.02 


12. 


,99 


R12B5X60 


220. 


,5 


2, 


.550 


0. 


,1770 


35.96 


2, 


.495 


9.983 


29.23 


17. 


,35 


R12B4X60 


223. 


,7 


3 


.011 


0. 


,01702 


43.08 


0, 


.2435 


7.955 


31.27 


11. 


,05 


R12B1X60 


217. 


,6 


1 


.933 


0. 


,2232 


27.24 


3 


.145 


12.26 


38.46 


15. 


,85 


R12B12X30 


230. 


,6 


2, 


.208 


0. 


,3719 


30.72 


5, 


,174 


17.36 


43.41 


14. 


,15 


R12B12X15 


209. 


,0 


1, 


.712 


0. 


,3392 


23.83 


4, 


,722 


16.61 


34.29 


28. 


,97 


R12B12X7 


231. 


,8 


1, 


.864 


0. 


,2876 


23.79 


3, 


.669 


110.0 


37.51 


34. 


,89 



